//
// Created by daiyan on 2021/6/17.
//
#include "Subdivide.h"
#include <OpenMesh/Core/Utils/PropertyManager.hh>

void loop_subdivide(OpenMesh::TriMesh& mesh)
{
    mesh.request_edge_status();
    mesh.request_face_status();
    mesh.request_vertex_status();

    using TriMesh      = OpenMesh::TriMesh;
    using Point        = TriMesh::Point;
    using EdgeHandle   = TriMesh::EdgeHandle;
    using VertexHandle = TriMesh::VertexHandle;

    auto newPosProp  = OpenMesh::makeTemporaryProperty<VertexHandle, Point>(mesh);
    auto newVertProp = OpenMesh::makeTemporaryProperty<EdgeHandle, VertexHandle>(mesh);

    // add edge middle point
    for (auto e : mesh.edges())
    {
        auto vHandle   = mesh.new_vertex();
        newVertProp[e] = vHandle;

        if (e.is_boundary())
        {
            // v0
            {
                auto        h0 = e.h0();
                auto        v0 = h0.from();
                const auto& p0 = mesh.point(v0);
                const auto& p1 = mesh.point(h0.to());
                const auto& p2 = mesh.point(h0.prev().from());
                newPosProp[v0] = 0.75 * p0 + 0.125 * (p1 + p2);
            }

            // vHandle
            {
                const auto& p0      = mesh.point(e.v0());
                const auto& p1      = mesh.point(e.v1());
                newPosProp[vHandle] = (p0 + p1) * 0.5;
            }
        }
        else
        {
            //
            {
                auto   v     = e.v0();
                auto   verts = v.vertices().to_vector();
                auto   n     = verts.size();
                double beta  = (0.625 - pow(0.375 + 0.25 * cos(2 * fPi / n), 2)) / n;

                Point p;
                p.vectorize(0);
                for (const auto& vi : verts)
                {
                    p += mesh.point(vi);
                }
                p *= beta;
                p += (1 - n * beta) * mesh.point(v);

                newPosProp[v] = p;
            }

            // vHandle
            {
                auto        h0      = e.h0();
                auto        h1      = e.h1();
                const auto& p0      = mesh.point(h0.from());
                const auto& p1      = mesh.point(h1.from());
                const auto& p2      = mesh.point(h0.next().to());
                const auto& p3      = mesh.point(h1.next().to());
                newPosProp[vHandle] = 0.375 * (p0 + p1) + 0.125 * (p2 + p3);
            }
        }
    }

    // collect new faces and delete old faces
    std::vector<std::vector<VertexHandle>> newFaces;
    for (auto f : mesh.faces())
    {
        std::vector<VertexHandle> V;
        for (auto he : f.halfedges())
        {
            V.push_back(he.from());
            V.push_back(newVertProp[he.edge()]);
        }

        /*
             0
             /\
           1/__\5
           /\  /\
         2/__\/__\4
              3
        */

        newFaces.push_back({V[0], V[1], V[5]});
        newFaces.push_back({V[1], V[2], V[3]});
        newFaces.push_back({V[1], V[3], V[5]});
        newFaces.push_back({V[5], V[3], V[4]});
        mesh.delete_face(f, false);
    }

    // assign new faces
    for (const auto& nf : newFaces)
    {
        mesh.add_face(nf);
    }

    // replace by properties
    for (auto i : mesh.vertices())
    {
        mesh.set_point(i, newPosProp[i]);
    }

    mesh.garbage_collection();
}

void square_root_3_subdivide(OpenMesh::TriMesh& mesh)
{
    throw Exception("not implement");
}
